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SUMMARY/OVERVIEW 

Practical  fuels,  such  as  JP8,  derived  from  the  refinery  processing  of  crude  petroleum,  are 
chemically  complex,  often  containing  thousands  of  compounds.  Because  these  fuels  need  to  meet 
broadly  defined  specifications,  their  compositions  vary  not  only  with  refinery  and  crude  oil 
source,  but  also  with  season  and  year  of  production.  As  a  result,  it  is  difficult  to  control  the 
consistency  in  fuel  composition  required  for  the  purpose  of  research.  When  dealing  with  these 
complex  fuels,  surrogates,  mixtures  of  a  limited  number  of  hydrocarbons  with  a  well-defined  and 
reproducible  composition,  can  be  used  in  place  of  the  real  fuel  both  for  experimental  and 
computational  applications. 

The  work  performed  by  our  research  group  over  the  last  years  aimed  to  develop  an  integrated 
and  innovative  computational  approach  for  predicting  reaction  mechanisms  for  JP8  surrogate 
components.  The  modeling  effort  has  integrated  various  computational  tools  to  build  chemical 
kinetic  pathways,  starting  from  the  structure  of  the  proposed  fuel  components  and  ending  with  a 
list  of  reactions  pathways,  rate  constants,  thermodynamic  and  transport  data. 

Studies  were  conducted  in  several  research  areas,  all  of  which  were  directed  at  meeting  the  overall 
project  objectives.  The  chemistry  of  H  abstraction  reactions  from  various  hydrocarbons  was 
obtained  using  ab  initio  methods.  Reaction  class  theory  was  applied  to  extend  the  knowledge 
obtained  with  quantum  chemistry  methods  to  large  classes  of  reactions.  Transport  properties  of 
n-alkane  were  determined  using  molecular  dynamics  simulations  in  conjunction  with  the  Green- 
Kubo  theory.  The  results  have  been  integrated  in  the  first  generation  of  a  jet  fuel  mechanism. 


TECHNICAL  DISCUSSION 

The  innovative  aspect  of  the  modeling  effort  carried  out  over  the  last  years  has  been  the 
integration  of  computational  tools  to  build  reaction  mechanisms  for  real  fuel  mechanisms, 
including  chemical  pathways,  rate  constants,  and  transport  data.  The  results  described  below  are 
organized  in  wo  main  areas:  kinetics  and  transport  properties  of  JP8  surrogate  components. 
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1.  H  abstraction  reactions:  ab  initio  simulations  and  Transition  State  Theory 

Hydrogen  abstraction  reactions  are  an  important  part  of  the  kinetic  mechanisms  available  in  the 
literature  to  study  combustion  systems.  We,  therefore,  focused  our  effort  on  obtaining  accurate 
rate  constants  for  various  hydrocarbons  using  multiple  radicals,  such  as  H,  OH,  CH3,  and  O,  as 
incoming  species  for  the  abstraction.  The  discussion  below  summarizes  the  results  obtained  for 
alkene  compounds. 

The  hydrogen  abstraction  reaction  between  a  hydroxyl  radical  (OH)  and  an  alkene  (C=C)  to  form 
a  water  molecule  (H20)  and  an  alkenyl  radical  (C=C#)  is  known  to  be  an  important  reaction  class 
in  combustion  processes  of  hydrocarbon  fuel,  especially  in  the  high  temperature  regime.  The 
hydrogen  abstraction  reaction  between  C2H4  and  OH  to  form  C2H3  and  H20  has  attracted  a 
number  of  extensive  experimental,  as  well  as  theoretical,  investigations.  There  are  more  than  15 
entries  for  rate  constant  studies  for  this  reaction  in  the  NIST  chemical  kinetics  database.  For 
reactions  involving  alkenes  larger  than  C2H4,  even  fewer  data  are  available  due  to  the  involvement 
of  other  kinds  of  reactions  such  as  the  addition  of  OH  to  the  double  bond  and  the  hydrogen 
abstraction  at  different  carbon  sites,  e.g.  saturated  carbon  sites  (sp  hybridization). 

All  the  electronic  structure  calculations  were  carried  out  using  the  GAUSSIAN  03  program. 
Hybrid  non-local  Density  Functional  Theory  (DFTQ,  particularly  Becke’s  half-and-half  (BH&H) 
non-local  exchange  and  Lee-Yang-Parr  (LYP)  non-local  correlation  functionals,  has  been  found 
to  be  sufficiently  accurate  for  predicting  the  transition  state  properties,  e.g.  barrier  height  and 
vibrational  frequency,  for  hydrogen  abstraction  reactions  by  a  radical.  Geometries  of  reactants, 
transition  states,  and  products  were  optimized  at  BH&HLYP  level  of  theory  with  the  Dunning’s 
correlation-consistent  polarized  valence  triple-zeta  basis  set,  denoted  as  cc-pVTZ,  that  is 
sufficient  to  capture  the  physical  change  along  the  reaction  coordinate  for  this  type  of  reaction. 
Frequencies  of  the  stationary  points  were  also  calculated  at  the  same  level  of  theory.  This 
information  was  used  to  derive  the  RC-TST  factors  described  in  the  following  section.  Thermal 
rate  constants  were  computed  for  the  temperature  range  of  300-3000  K.  The  correction  for  the 
hindered  rotors,  e.g.  reduced  moment  of  inertia  and  the  periodic  potential,  are  then  obtained 
using  the  information  on  the  rate  constants. 

Figure  1  reports  the  rate  constants  available  in  the  literature  for  this  reaction  obtained  by 
experiments  and  simulations.  Tully  et  al.  1  used  the  laser  photolysis/laser-induced  fluorescence 
technique  under  slow- flow  conditions  to  measure  the  rate  constants  in  the  temperature  range  of 
650-901  K.  Baulch  et  al.  presented  critically  evaluated  kinetic  data  for  use  in  computer 
combustion  modeling.  The  suggested  rate  data  with  an  uncertainty  factor  of  3  in  the  temperature 
range  of  650-1500  K  follow  closely  the  studies  of  Tully  et  al.  x  Westbrook  et  al.  3  suggested  rate 
constant  data  for  the  reference  reaction  in  the  temperature  range  1003-1253  K.  Using  visible-UV 
absorption  technique  together  with  TST  model,  Liu  et  al.  4  reported  the  rates  of  the  hydrogen 
abstraction  in  the  temperature  range  of  723  -  1170  K.  With  the  energetic  data  obtained  at  the 
QCISD(T)/6-311G(2df,p)//B3LYP/6-31G(d,p)  together  with  the  gradient  and  Hessian 
information  at  the  B3LYP/6-31G(d,p),  Liu  et  al.5  reported  thermal  rate  constants  in  the 
temperature  range  of  200-5000  K  using  the  canonical  variational  transition  state  theory  (CVT) 
and  the  small-curvature  tunneling  correction  (SCT).  Recently,  Senosiain  et  al.6  also  carried  out  an 
analysis  of  the  0H+C2H4“^H20+C2H3  reaction  using  the  CVT  method  with  the  molecular- 
property  data  obtained  at  the  RQCISD(T)/cc-pVooZ//UB3LYP/6-31 1  +  +G(d,p)  level  of  theory 
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to  suggest  rate  constants  for  this  reaction.  The  two  rate  constants  obtained  with  the  CVT  theory 
are  similar  (within  a  deviation  factor  of  4)  in  the  high  temperature  regime  but  they  differ  by  an 
order  of  magnitude  at  lower  temperature  (<  500  K). 


0  0  0  5  1  0  1  5  2  0 
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F  gure  1  Arrhenius  plots  of  rate  constants  for  the  reaction 
OH*C;H4->  H>0*C2Hj.  The  error  bars  tor  these  rate  constants  are 
also  included 


The  rate  constant  derived  in  our  study  for  the  reaction  OH  +  C2H4  -)  H20  +  C2H3  has  the 
following  expression: 

[cmVnolecule  V1] 

In  a  similar  way,  we  determined  rate  constants  for  H  abstraction  reactions  from  other  key 
reactions,  including  aromatics  containing  5  and  6-member  rings. 

2.  Reaction  Class  Theory  for  classes  of  hydrocarbons 

Detailed  chemical  kinetic  descriptions  of  real  fuel  combustion  require  the  tracking  of  hundreds  of 
chemical  species  and  thousands  of  reaction  steps.  It  is  impracticable  to  obtain  the  correct  kinetic 
data  for  such  a  large  number  of  reactions  by  experiments  or  explicit  rate-constant  calculations 
even  using  the  simple  Transition  State  Theory  method.  To  this  end,  studies  were  conducted  bv 
our  research  group  using  the  Reaction  Class  Transition  State  Theory  for  esumating  rate  constants 
of  a  large  number  of  reactions.  As  an  example,  below  we  report  the  results  obtained  for  the 
OH+alkene  reaction  class.  The  main  task  was  to  find  correlation  expressions  between  rate 
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constants  of  the  reference  reaction  (C2H4  4-  OH  =  C2H3  4-  H20)  and  those  of  other  reactions  in 
the  class  (CnH2n  4-  OH)  from  explicit  direct  ab  initio  dynamics  calculations  of  the  rate  constants  of 
a  representative  set  in  the  class  itself. 

To  compute  the  Reaction  Class  Transition  State  Theory  (RC-TST)  parameters  for  the 
OH+alkene  class,  15  reactions  are  considered  as  a  representative  set.  These  reactions  are  given 
below: 
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where  tram  and  as  denote  trans-  and  cis-  configurations  for  the  carbon  chain.  Note  that  this  set 
does  not  include  reactions  with  resonance  systems,  e.g.  1,3-butadiene,  as  well  as  aromatic 
systems,  e.g.  benzene. 

The  basic  idea  of  the  RC-TST  technique  is  that  reactions  belonging  to  a  specific  class  have  the 
same  reactive  moiety;  thus  the  difference  between  the  rate  constants  of  any  two  reactions  in  the 
class  is  mainly  due  to  the  difference  in  the  interactions  between  the  reactive  moiety  and  their 
different  substituents.  Within  the  RC-TST  framework,  the  rate  constants  of  an  arbitrary  reaction 
(denoted  as  kj  are  proportional  to  those  of  a  reference  reaction,  kn  (usually  the  smallest  reaction 
in  the  class,  which  is  referred  to  as  the  principal  reaction)  in  the  same  class  by  a  temperature 
dependent  function  f(T) : 

ka  (T)  =  f(T)  x  lcr  (T)  (1) 

The  rate  constants  for  the  reference  reaction  are  often  known  experimentally  or  can  be  calculated 
accurately  from  first-principles.  The  key  idea  of  the  RC-TST  method  is  to  factor  f(T)  into 
different  components  under  the  TST  framework: 

fcn-fa  x/«  xfgxfv  (2) 

where  fa ,  fK ,  fQ  and  fv  are  the  symmetry  number,  tunneling,  partition  function  and  potential 

energy  factors,  respectively.  These  terms  are  simply  the  ratios  of  the  corresponding  components 
in  the  TST  expression  for  the  arbitrary  and  reference  reactions: 
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where  O  is  the  reaction  symmetry  number;  k(T)  is  the  transmission  coefficient  accounting  for 
the  quantum  mechanical  tunneling  effects;  Q *  and  O*  are  the  total  partition  functions  (per  unit 
volume)  of  the  transition  state  and  reactants;  AV *  is  the  classical  reaction  barrier  height;  T is  the 

temperature  in  Kelvin;  kB  and  h  are  the  Boltzmann  and  Planck  constants,  respectively.  The 
potential  energy  factor  can  be  calculated  using  the  reaction  barrier  heights  of  the  arbitrary 
reaction  and  the  reference  reaction.  The  classical  reaction  barrier  height  AV *  for  the  arbitrary 

reaction  can  be  obtained  using  the  Linear  Energy  Relationship  (LER),  similar  to  the  well-known 
Evans-Polanyi  linear  free  energy  relationship,  '  between  classical  barrier  heights  and  reaction 
energies  in  a  given  class  without  having  to  calculate  them  explicidy.  Alternatively,  the  barrier 
height  for  the  arbitrary  reaction  can  be  obtained  from  the  Barrier  Height  Group  (BHG)  approach 
where  reactions  in  a  sub-class  can  be  reasonably  assumed  to  have  the  same  barrier  height. 


The  potential  energy  factor  can  be  calculated  using  Eq.  (6),  where  AV*  and  AV*  are  the  barrier 

heights  of  the  arbitrary  and  reference  reactions,  respectively.  We  have  also  shown  that  within  a 
given  class  there  is  a  linear  energy  relationship  (LER)  between  the  barrier  height  and  the  reaction 
energy,  similar  to  the  well-known  Evans-Polanyi  linear  free  energy  relationship.  Thus,  accurate 
barrier  heights  can  be  predicted  from  only  the  reaction  energies.  The  barrier  heights  for  reactions 
RrR15  can  also  be  grouped  into  two  classes:  terminal  carbon  sites  of  the  double  bond  (class  1)  and 
non-terminal  carbon  sites  (class  2).  This  can  be  referred  to  as  the  barrier  height  grouping  (BHG). 
The  AMI  semi-empirical  method  was  employed  to  calculate  the  reaction  energies.  AMI  and 
BH&HLYP/cc-pVTZ  reaction  energies  were  then  used  to  derive  the  LER’s  between  the  barrier 
heights  and  reaction  energies. 

The  observed  LERs  plotted  against  the  reaction  energies  calculated  at  BH&HLYP/cc-pVTZ  and 
AMI  levels  are  shown  in  Figs.  2a  and  2b,  respectively.  The  substitute  of  an  alkyl  group  will 
stabilize  the  radical  products  thus  lowering  the  barrier  heights.  For  this  reason  the  reactions  at 
the  non-terminal  carbon  of  the  double  bond  (class  2)  have  barrier  heights  of  about  2.0  keal/mol 
lower  than  those  at  the  terminal  sites. 

The  linear  fits  for  these  reactions  were  obtained  using  the  least  square  fitting  method  and  have 
the  following  expressions: 
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A V"  =0.4892  xAEbh&hlyp  +  10.772  (keal/mol)  (8a) 

AV  =  0.4238  xAfi"*"1  +16.572  (keal/mol)  (8b) 

Based  on  the  BHG  results,  we  assigned  the  values  of  10.50  and  8.27  keal/mol  to  the  energy 
barriers  of  the  terminal  and  non- terminal  carbon  sites  of  the  double  bond,  respectively.  It  is 
interesting  to  note  that  the  averaged  deviations  of  reaction  barrier  heights  estimated  from  the 

BHG  (0.15  keal/mol)  is  smaller  than  that 
of  the  LER,  while  the  maximum  deviation 
(0.46  keal/mol)  is  higher.  Therefore,  this 
approach  can  be  used  to  estimate  the 
relative  barrier  height  quickly  with  an 
acceptable  confidence.  The  key  advantage 
of  this  approach  is  that  it  does  not  require 
any  additional  information  to  estimate 
rate  constants. 

In  conclusion,  the  barrier  height  of  any 
reaction  in  the  OFI+alkene  reaction  class 
can  be  obtained  by  using  either  the  LER 
or  BHG  approach.  The  estimated  barrier 
height  is  then  used  to  calculate  the 
potential  energy  factor  using  Eq.  (6).  The 
performance  of  both  approaches  is 
discussed  in  the  error  analyses  below. 

The  tunneling  factor  JK  is  the  ratio  of  the 
transmission  coefficient  of  reaction  Ra  to 
that  of  reaction  R^.  Due  to  the 
cancellation  of  errors  in  evaluating  the 
tunneling  factors,  we  have  shown  that  the 
factor  fK  can  be  reasonably  estimated 
using  the  one-dimension  Eckart  method. 
In  the  Eckart  formulation,  the  imaginary 
frequency  and  the  barrier  height  are  used  to  calculate  the  tunneling  probability  for  a  reaction. 
Because  the  barrier  heights  are  grouped  into  2  classes,  namely  terminal  and  non-terminal  sites  of 
the  double-bond  carbon  and  the  imaginary  frequencies  for  reactions  at  the  same  class  are  very 
similar,  the  values  of  the  tunneling  coefficients  for  reactions  in  the  same  class  are  expected  to  be 
similar.  Therefore,  the  average  value  for  the  tunneling  factors  can  be  used  for  the  whole  group. 
Simple  expressions  for  the  two  tunneling  factors  for  terminal  and  non-terminal  carbon  sites  of 
the  double  bond  are  obtained  by  fitting  to  the  average  calculated  values  and  are  given  below: 

fK  =  0.999  -  83 .42  x  exp[-0.5 1  x  T° 42 ]  for  terminal  carbon  sites  (9a) 

fK  =  0.978  -  7.55  x  exp[-0.056x  T{  66]  for  non-terminal  carbon  sites  (9b) 

The  correlation  coefficients  for  these  fits  are  larger  than  0.999.  The  same  tunneling  factor 
expression  can  be  reasonably  assigned  to  different  reactions  in  the  same  class  with  the  largest 


Figure  2:  Linear  energy  relationship  plots  of  the  barrier 
heights  versus  the  reaction  energies,  computed  at  (a) 
BH&HLYP/cc-pVTZ  and  (b)  AMI  levels  of  theory. 
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unsigned  deviation  of  0.08  for  R13  and  the  largest  percentage  deviation  of  17.2%.  The  mean 
unsigned  deviation  is  7%,  compared  to  the  direct  Eckart  calculation  using  reacuon  information 
from  BH&HLYP/cc-pVTZ  level  of  theory.  At  higher  temperatures,  tunneling  contributions  to 
the  rate  constants  decrease  and  thus,  as  expected,  the  differences  between  the  approximated 
values  and  the  explicidy  calculated  ones  also  decrease;  for  example,  the  maximum  error  for  all 
reacuons  is  less  than  2%  at  500  K. 

The  partition  factor  includes  translational,  rotational,  internally-rotational,  vibradonal,  and 
electronic  components.  The  partidon  funcdon  factor  fg  mainly  originates  from  the  differences 

in  the  coupling  between  the  substituents  and  the  reacdve  moiety.  Its  temperature  dependence 
arises  from  the  vibradonal  and  internally-rotadonal  components  only. 

The  averaged  values  of  partition  function  factors  for  reacdons  R2-R15  are  fitted  into  the  following 
analydcal  expression: 

fQ  =0.71-2.08xexp[-0.18x7’°45]  (10) 

The  coupling  between  subtituents  with  the  reacdve  moiety  is  believed  to  account  for  the  partition 
funcdon  factors  having  values  of  around  0.7.  The  total  coupling  effect  is  contributed  from  those 
of  the  translational,  rotational  and  vibradonal  partition  function  factors.  Each  reaction  class  has  a 
specific  coupling  effect  mainly  due  to  the  specific/unique  reactive  moiety.  If  there  is  no  coupling 
effect,  the  values  of  the  partidon  funcdon  factors  would  be  expected  to  be  very  close  to  unity. 

For  this  reaction  class,  the  rotation  of  the  alkyl  group  (CH3)  along  the  C — C  bond  at  some 
reactants,  transition  states  and  products  as  well  as  the  rotation  of  the  hydroxyl  (OH)  group  along 
the  C — O  axis  at  all  transition  states  need  to  be  treated  as  hindered  rotations.  The  hindered 
rotation  treatment  lowers  the  total  rate  constants  with  the  temperature  increase.  Note  that  the 
principal  reaction  R t  does  not  have  the  internal  rotation  of  the  CH3  group.  The  averaged  values, 
as  applied  to  the  whole  class,  are  fitted  into  an  analytical  expression  as  given  below: 

fHR  =  1.0 1-0.72  xexp[- 1332  xT'° 99  ]  (11) 

So  far  we  have  established  the  necessary  parameters  -  namely  the  potential  energy,  the  symmetry 
number,  the  tunneling  and  the  partition  function  factors  -  for  application  of  the  RC-TST  theory 
to  predict  rate  constants  for  reactions  in  the  OH+Alkene  class. 

The  procedure  for  calculating  the  rate  constants  of  an  arbitrary  reaction  in  this  class  is:  a) 
calculate  the  potential  energy  factor  using  Eq.  (6)  with  the  AV*  value  of  11.07  kcal/mol.  The 
reaction  barrier  height  can  be  obtained  using  the  LER  approach  by  employing  Eq.  (8a)  for 
BH&HLYP/cc-pVTZ  reaction  energies  or  (8b)  for  AMI  reaction  energies  or  by  the  BHG 
approach;  (b)  calculate  the  symmetry  number  factor  from  Eq.  (3)  or  see  Table  2;  (c)  compute  the 
tunneling  factor  using  Eq.  (9a)  and  (9b)  for  terminal  and  non-terminal  carbon  sites,  respectively; 
(d)  evaluate  the  partition  function  factor  using  Eq.  (10)  with  the  hindered  roration  treatment 
correction  using  Eq.  (1 1);  and  (e)  the  rate  constants  of  the  arbitrary  reaction  can  be  calculated  by 
taking  the  product  of  the  reference  reaction  rate  constants  given  by  Eq.  (7)  with  the  reaction  class 
factors. 
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Table  2  summarizes  the  RC-TST  parameters  for  this  reaction  class.  If  the  BHG  barrier  heights 
and  the  average  values  for  the  factors  are  used,  the  rate  constants  are  denoted  by  RC-TST/BHG. 
The  RC-TST /BHG  rate  constants  for  any  reaction  belonging  to  this  reaction  class  are: 


k(T)  *  7.51x10 


-24 


k(T)  =  1.00x10” 


rr>3.59 

xT  x  exp 

'185' 

T 

ml  19 

'763' 

xT  x  exp 

T 

for  terminal  carbon  sites 


for  non-terminal  carbon  sites 


(12a) 

(12b) 


Because  the  terminal  carbon  sites  have  two  hydrogen  atoms  which  can  be  reasonably 
considered  equivalent  in  some  cases,  and  the  non-terminal  sites  only  have  one  hydrogen  atom, 
the  symmetry  factors  of  2  and  1  are  also  included  in  the  rate  constant  expressions. 


TABLE  2:  Parameters  and  Formulations  of  the  RC-TST  Method  for  the  OH  +  Alkene 
H20  +  Alkenyl  Reaction  Class  (OH+  C2H4  is  the  reference  reaction) 


k(T)-f„x  fK  (T)  x  f0  (T)  X  fHR(T)  X  fv  (T)  x  k,(T) ;  fv  (. T )  =  exp 


-(AF*-AF*) 


kBT 


T  is  in  Kelvin;  AV*and  AE  are  in  kcal/mol;  Zero-point  energy  correction  is  not  included 
fa  Calculated  explicitly  from  the  symmetry  of  reactions  (see  Table  2) 


fK  =  0.999  -83.42 xexp[-0.51  xT"4  ]  for  terminal  carbon  sites. 


/AT) 

fK  =  0.978  -  7.55  x  exp[— 0.056  x  T° 66 ]  for  non-terminal  carbon  sites. 


/q(T)  fo  =  0.71  -2.08  xexpI-O.lSxT0  45  ] 


IhAT) 

fHR  =  1 .0 1  -  0.72  X  exp[- 1 332  X  r° 99  ] 

AFJ 

A Vx  -0.4892  xAEbh&hlyp  +  10.772 

LER 

AF*  =  0.4238  xAf^1  +16.572 

AFr*  =-11.07  kcal/mol 

kr(T) 

(Eq.7) 

kr  (T)  =  2.18x1 0"25  x  T4  20  x  exp 

14331 

T 

,  [cm^molecule  ls  l] 
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To  illustrate  the  theory  we  selected  two 
reactions  R2  and  R3  whose  rate  constants 
are  available  in  the  literature. 

Figs.  3a~b  show  the  predicted  rate 
constants  of  reaction  R2  and  reaction  R3 
using  the  RC-TST  method.  Since  the 
barrier  heights  obtained  from  either 
BH&HLYP/cc-pVTZ  or  AMI  energies 
are  similar,  we  can  expect  their  rate 
constants  to  be  similar. 

The  rate  constants  estimated  from  the 
RC-TST/LER  and  RC-TST/BHG 

approaches  are  comparable  for  reactions 
R2  and  R3  due  to  the  similar  predicted 
values  of  the  barrier  heights,  e.g.  10.47 
and  10.50  kcal/mol  for  R2  from  LER  and 
BHG,  respectively.  The  difference  in  rate 
constants  might  be  larger  for  other 
reactions. 

Compared  to  the  “RC-TST  exact”  values, 
the  excellent  performance  of  the  RC-TST 
for  these  two  reactions  can  be  seen  in 
Figure  3. 
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Figure  3:  Arrhenius  plots  of  the  calculated  rate  constants 
using  the  RC-TST  method  for  two  representative  hydrogen 
abstraction  reactions  along  with  the  available  values;  (a) 
OII+C31I6  at  the  terminal  carbon,  and  (b)  OII+C3116  at  the 
non-terminal  carbon. 


3.  Diffusion  Coefficients  of  Fuel  Surrogate  Components  using  Molecular  Dynamics 
Simulations 

In  solving  chemically  reacting-flow  problems,  however,  chemical  production  and  destruction  is 
often  balanced  by  transport  due  to  convection,  diffusion,  or  conduction  and,  hence,  accurate 
transport  coefficients  are  of  great  importance  in  the  development  of  reaction  mechanisms. 
Therefore,  in  addition  to  kinetic  studies,  we  have  initiated  an  analysis  of  the  transport  properties 
of  hydrocarbons  using  Molecular  Dynamics  simulations.  Specifically,  we  have  looked  at  the 
binary  diffusion  coefficients  of  n-alkane  molecules. 

The  Chapman  Enskog  (CE)  relationship  is  the  theoretical  equation  currently  used  to  compute  the 
diffusion  coefficients  of  different  systems.10, 11  Combined  with  the  Lennard-Jones  intermolecular 
potential  function,  the  binary  diffusion  coefficient  can  be  expressed  as: 

D  3  V2jr(^n3/m,2 

12  16  Pjio?2Q 
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where  kB  is  the  Boltzmann  constant  and  ml2  is  the  reduced  mass  of  the  pair  components.  This 
expression  contains  a  dependence  on  the  temperature,  T,  pressure,  P,  the  average  collision 
diameter  between  two  species,  On ,  and  the  collision  integral,  Q.  The  collision  integral  term  is 
related  to  the  interaction  energy  of  two  molecules  and  the  scattering  phenomena  when  the 


collision  u  The  value  of  the  collision  integral  depends  on  the  reduced  temperature,  T 


where  e]2  is  the  well  depth  of  the  Lennard-Jones  energy  potential. 

This  approach  presents  many  limitations.  The  CE  equation  considers  the  gases  as  hard  spheres 
and  requires  accurate  values  of  the  collision  diameter  and  potential  energy  well  depth.  These 
parameters  are  usually  not  available  and  in  many  cases  they  are  assumed  based  on  similarity. 


The  CE  equation  does  not  give  reliable  transport  data  under  high  density  or  high  temperature 
conditions  especially  for  binary  diffusion  case.  The  Chapman-Enskog  equation  is  valid  only  for 

p  <0.1  where 

P  =  P°\ 

p  is  the  number  density  of  the  system  and  Oe  is  the  effective  hard  sphere  diameter  which 
depends  on  temperature. 13 

Finally,  the  CE  equation  does  not  take  into  account  the  effect  of  molecular  geometry  variations 
on  the  mutual  diffusion  coefficient  of  isomers,  an  important  factor  especially  in  high  density  or 
high  mobility  condidons. 


The  Green-Kubo  relations  provide  the  exact  mathematical  expression  for  transport  coefficients 
in  terms  of  microscopic  fluctuations  of  the  system  of  interest.14  For  the  self  diffusion  coefficient, 
the  Green  Kubo  is  the  time  integral  of  the  velocity  auto-correlation  function,  which  represents 
the  correlation  between  molecules  at  different  times1" 


A  = 


3  Nt 


dt 


where  N]  is  the  total  number  of  particles  of  species  1,  Uj  is  the  velocity  component  of  particle  i 
of  species  1,  and  the  angular  brackets  denote  the  ensemble  average  of  velocity  auto  correlation 
functions.  The  binary  diffusion  coefficient  can  be  expressed  as  the  combination  of  the  time 
integral  of  the  velocity  auto-  and  cross-correlation  function,  which  is  the  correlation  between  the 
different  molecules  at  different  times. 


Dn-Q 


x2D]  +  x]D2  +  xxx2  f  +  2-^ 


*1  X2  XIX2 


l  He. 


where  Q  is  the  thermodynamic  factor  determined  from  the  integral  of  the  radial  distribution 
functions  with  respect  to  the  ensemble,16  x{  and  x2  are  the  mole  fractions  of  the  species  1  and  2; 
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Dj  and  D2  arc  the  self  diffusion  coefficients  of  the  species  1  and  2  in  the  mixture;  fap  is  the 
ensemble  average  of  the  time  integral  of  the  velocity  cross-correlation  function 

Molecular  Dynamics  simulations  have  been  used  to  denve  the  microscopic  information  in  terms 
of  the  time  integral  of  the  velocity  cross-correlation  components  of  the  species  in  various 
temperature  and  pressure  regimes.  As  initial  ensemble,  we  studied  n-alkane  compounds  with  the 
goal  to  determine  the  accuracy  of  C-E  in  dealing  with  long  molecules. 
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Figure  4:  Binary  mass  diffusivities  of  CH4 
(upper  panel)  and  C16H34  in  nitrogen  as 
function  of  temperature.  CE:  Chapman-Enskog; 
MD:  Molecular  Dynamics. 


Figure  4  shows  the  comparison  between  the  binary  mass  diffusivities  of  methane  and  hexadecane 
in  nitrogen  obtained  using  the  C-E  equation  and  MD  simulations.  While  the  agreement  with  the 
smallest  molecule  (almost  a  spherical  one)  is  very  good,  the  comparison  for  the 
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hexadecane/nitrogen  mixture  (long  chain  structure)  shows  a  large  deviation  reaching  17.3  %  at 
high  temperatures. 

The  table  below  shows  the  relative  deviations  between  the  diffusion  coefficients  obtained  using 
C-E  and  MD  for  the  series  of  alkane  Cl -Cl  6.  For  mixtures  composed  of  hydrocarbons  in  the 
range  C1-C6  the  deviation  is  less  than  5%.  The  MD  results  consistently  show  lower  diffusivities  than 
those  obtained  using  the  C-E  equation  due  to  the  anisotropic  force  field  interactions. 


C2H/S2 

8.627X1  (rT1-714 

c3h8/n2 

6.479x1 0’V’725 

c4h10/n2 

5.189x10'6T1735 

c*h12/n2 

4.361  xlO^T1 743 

c6h14/n2 

3.449x1 0'V'750 

c7h16/n2 

3.059x104’T1,757 

C8Hig/N2 

2.706x1  O^T1'764 

Ci0H22/N2 

2.124x10‘‘t1'775 

c12h26/n2 

1.740xl0‘"T1-784 

C14H30/N2 

l^xlO"4'!'1,7*2 

Ci6H34/N2 

1.309xlO‘‘T1796 

Table  3.  Binary  mass  diffusion  coefficients 
of  n-alkanes  in  nitrogen  as  a  function  of 
temperature  [K]  at  1  atm. 


4.  Integration  of  the  various  projects 

Results  from  the  projects  listed  above  have  been  implemented  in  an  initial  version  of  a  kinetic 
mechanism  for  JP8  surrogate  components  (JetSurf).  The  model  is  being  developed  through  a 
multi-university  research  collaboration  involving  experimental  and  modeling  groups  at  the 
University  of  Southern  California,  Stanford  University,  NIST,  Drexel  University,  Princeton 
University  and  the  Imperial  College,  London. 
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